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Abstract 

We consider in this paper, a few important issues in non-equilibrium work fluctuations and 
their relations to equilibrium free energies. First we show that Jarzynski identity can be viewed 
as a cumulant expansion of work. For a switching process which is nearly quasistatic the work 
distribution is sharply peaked and Gaussian. We show analytically that dissipation given by average 
work minus reversible work Wr, decreases when the process becomes more and more quasistatic. 
Eventually, in the quasistatic reversible limit, the dissipation vanishes. However estimate of p - the 
probability of violation of the second law given by the integral of the tail of the work distribution 
from —00 to Wr, increases and takes a value of 0.5 in the quasistatic limit. We show this analytically 
employing Gaussian integrals given by error functions and Callen-Welton theorem that relates 
fluctuations to dissipation in process that is nearly quasistatic. Then we carry out Monte Carlo 
simulation of non-equilibrium processes in a liquid crystal system in the presence of an electric 
field and present results on reversible work, dissipation, probability of violation of the second law 
and distribution of work 
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I. INTRODUCTION 



Non equilibrium work fluctuations and equilibrium free energy differences are related to 



each other through a remarkable identity proposed and proved by Jarzynsk 



2] . Ever since 



several studies have been carried out to establish the validity and utility of the identity for 



various non-equilibrium evolutions - stochastic, see e.g. 



BE 



4| as well as deterministic, see 



e.g. 

We show in this paper that Jarzynski identity can be viewed as a cumulant expansion. 
When we retain only the first cumulant in the expansion, we get the well known thermody- 
namic identity equating the work done on the system to free energy change in a quasistatic 
reversibleprocess. With only the first two cumulants in the expansion, we get Callen-Welton 
theorem [6j|. When we include the third and higher order cumulants we get Jarzynski iden- 
tity. We derive an analytical expression for probability of violation of the second law, under 
Gaussian approximation to work fluctuations. We then carry out Monte Carlo simulation 
of a liquid crystal system in the presence of a varying electric field and report numerical 
results on free energy changes, dissipation, probability of violation of the second law and 
distribution of work. 

The paper is organized as follows. We start with a brief introduction to a few relevant 
basic issues of heat and work in a quasistatic reversible process. Then we consider heat and 
work in the context of irreversible processes. This is followed by a discussion on the relation 
between non-equilibrium work fluctuations and equilibrium free energy differences. We show 
analytically that the probability of violation of the second law increases with increase of time 
duration of the switching experiment and in the asymptotic limit of a quasistatic process 
it goes to one-half. However dissipation defined as average work minus reversible work 
becomes smaller when the duration of time of the switching process increases and in the 
quasistatic limit it goes to zero. We demonstrate the usefulness of Jarzynski identity 
on a system of liquid crystals in the presence of external electric field. We consider a lattice 
model with each lattice site holding a headless spin. We calculate the work done when 
the external electric field is switched from an initial value to a final value following a well 
specified experimental protocol. We carry Monte Carlo simulation of the non-equilibrium 
switching process and collect an ensemble of work values. From this ensemble we extract 
equilibrium properties of the liquid crystalline system employing both Jarzynski identity 



and Callen-Welton theorem |6j. Principal conclusions of the study are brought out briefly 
in the concluding section. 



II. HEAT AND WORK : REVERSIBLE PROCESSES 

Consider a closed system in equilibrium at temperature T. It draws a small quantity 
d Q of heat by a quasistatic reversible process at constant temperature. dQ is not a perfect 
differential. However, 

(1) 

is a perfect differential where S is entropy. Thus, during a quasistatic reversible process, 
entropy of the system increases by ctQ/T. 

To calculate the work done on the system, we start with the first law of thermodynamics, 
stated as, 

dW = dU - dQ, (2) 

where dW denotes the work done on the system and U, the internal energy. dW is not a 
perfect differential. 

The first law is about conservation of energy and is valid for all processes, quasistatic, 
non-quasistatic, reversible, irreversible or otherwise. However, if the process is quasistatic 
and reversible, then we can replace dQ by TdS, see Eq. (TjQ), and write, 

dW = dU - TdS. (3) 

If the process is also isothermal, we have 

dW = d(U -TS). (4) 

We identify the term U — TS as Legendre transform of the fundamental equation U = 
U(S, V) where we transform the variable S in favor of the 'slope' 



ns,v) = {§) (5) 



and the variable U in favor of the 'intercept' 

F(T, V) = U(S, V) - TS. (6) 



F(T, V) is called Helmholtz free energy or simply free energy. Thus the work done equals 
free energy change : 

dW = dF. (7) 

For purpose of distinguishing process variables (heat and work) from state variables (energy 
U, entropy S and free energy F) we have considered in the above an infinitesimal process. 
dQ and dW are not perfect differentials; dll, dS and dF are perfect differentials. However 
the relation between free energy and work, see Eq. ([7]), holds good for a quasistatic reversible 
process that takes the system from one equilibrium state, say A to another equilibrium state, 
say B. We have 

W = F(B) - F(A) (8) 
We call the above as reversible work and denote it by Wr 



III. HEAT AND WORK : IRREVERSIBLE PROCESSES 



To obtain a relation between change in free energy and work done when the process is 



not reversible, we start with the second law inequality, see js] for an elegant proof, given by 

dS > (9) 

The above implies 

dQ < T dS (10) 

for an irreversible process. Substitute this in Eq. (T5|) and get, 

dW > dU -TdS (11) 

which for an isothermal irreversible process reduces to 

dW > dF. (12) 

Consider a process that takes the system irreversibly from an equilibrium state A to another 
equilibrium state B. Let W denote the work done during the process and AF = F(B) — 
F(A). We have 

W > AF (13) 



The work done exceeds free energy change, if the process is not reversible. This is the best 
conclusion we can draw from purely thermodynamics considerations. 

A natural question that arises in this context is about the meaning of the statement 
W > AF. A given process can be realized either in an experiment or in a computer 
simulation. We recognize that work done, in general, would differ from one experiment 
or computer realization of the process to another. Let us consider a process in which a 
parameter of the thermodynamics system is switched from one value to another, as per a 
well defined experimental protocol. For example, we can change the volume of a system from 
one value to another uniformly over a fixed duration of time r. We call this a switching 
process. In general, different switching experiments, all carried out with the same protocol, 
will give rise to different values of W. Only when the switching is done quasistatically and 
reversibly does one get the same value of W in all experiments. Hence, in general, we have 
to deal with an ensemble of values of W and not just a single value. Let us denote this 
ensemble by Q = {Wi}. It is then quite possible that there can exist realizations for which 
W is less than dF, thereby violating the second law. This should not surprise us since even 
in the early days of statistical mechanics, Maxwell correctly recognized that the second 
law is of statistical origin and hence there is a non-zero probability of it being contravened. 
The demon he created to drive home this point, haunts us even today, see e.g. [1 Of ] . 

Let the ensemble Q of work values be formally described by the distribution p(W; r) 
where r denotes the duration of time over which the process takes place; r is the switching 
time. Typically we switch a macroscopic parameter, denoted by the symbol A from an initial 
value say Aj to a final value say Ap. This switching can be carried out in any way. The 
discussions below and the Jarzynski identity described in the next section hold good for any 
protocol of switching. However we consider the switching to take place at constant rate, 
over a pre specified duration of time, r. Then by considering different values of r we can 
describe different switching scenarios: if r is small, we have fast switching; if r is large we 
have slow switching; in the limit of r — > oo, we have a quasistatic reversible switching. Thus 



The probability of violation of the second law, denoted by the symbol p(r) is formally given 
by the following integral. 




< t < r. 



(14) 



T 




(15) 



5 



We note that for every thermodynamic variable there corresponds a random variable in 
statistical mechanics. The average of the random variable over a suitable ensemble gives the 
value of the thermodynamic variable. For example (E) equals the thermodynamic energy 
U, where the angular brackets denote average over a canonical ensemble. Hence, strictly we 
should state the second law as 

( W ) > dF . (16) 

In the above, the angular brackets denote an average over an ensemble of switching experi- 
ments, all carried out with the same protocol. In other words it is given by 



(W) = dWW p(W;r) (17) 

J — oo 

Stated thus, the second law can never be violated. We define 

W d = (W) - W R (18) 

as dissipation and state the second law as described in the next section 

W d > (19) 

for all processes. In the above, equality obtains when the process is quasistatic and reversible. 

If r is very large but not infinity, p(W; r) would be sharply peaked and Gaussian. Consider 
a process that takes a system from an equilibrium state A to another equilibrium state B. 
Let the process be nearly quasistatic. We have AF = F(B) — F(A) = Wr. The time taken 
for the process, r, is large. For such a process, dissipation is proportional to fluctuations, 



and from Callen-Welton theorem jf}] we get 



where 



is the first cumulant and 



AF = W R = Ci - ~ d (2 (20) 



Ci = (W) (21) 



( 2 = a 2 w = (W 2 ) - (W) 2 (22) 



is the second cumulant (or variance) of W. Thus we have Wd = $(2/2 from Eq. (TSUI) : 
dissipation (Wd) is proportional to fluctuation (£ 2 or crj^). 
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Measuring energy in units of fcgT = (3 1 , we get, from Eq. ( )20|) . 

- /? dF = {-13) Ci + - (-/3) 2 C2 • (23) 

The above relation is for a process which is nearly quasistatic. The distribution of work for 
such a switching protocol is Gaussian. For a Gaussian, the third and higher order cumulants 
are all identically zero. We contend that in a general process, which is not necessarily 
quasistatic or near quasistatic, there should be additional terms in the Right Hand Side 
(RHS) of the above equation involving third and higher cumulants. Once we include these 
higher cumulants we get Jarzynski identity, as shown in the next section. 

IV. NON-EQUILIBRIUM WORK AND EQUILIBRIUM FREE ENERGIES 

Including the terms involving the third and higher order cumulants in the RHS of Eq. 
(E3D, we get, 

-0dF= (-IS) Ci + C 2 + £ ^ Cn- (24) 

l - n=3 U - 

We recognize immediately the Right Hand Side (RHS) of the above equation as the cumulant 
expansion of W. 

Let x(P) denote the cumulant generating function. Thus we have, 

X ((3) = -(3dF (25) 

or equivalently 

exp(x) = exp(-/3dF). (26) 
The moment generating function of W is defined as 

0(/3) = / dW exp(-pW) p(W;t) 



= (exp(-/3W)). (27) 
The moment and cumulant generating functions are related to each other. We have 

X{P) = log W)\ ■ (28) 
We see immediately that Jarzynski identity given by, 

( exp(-pW) ) = exp(-pdF), (29) 

follows naturally from this. 



V. GAUSSIAN WORK DISTRIBUTION 



We shall derive an expression for the probability of violation of the second law formally 
given by Eq. (}15|) . For large r the work distribution p(W; r) would be a sharply peaked 
Gaussian with mean (W(r)) and variance a w (r) = (W /2 (r)) — (H^r)) 2 . Substituting the 
Gaussian in Eq. (ITS]) for p(W / ; r) we get, 



p(r) = - — / dW exp 1 



2> 



cr w (r)A/27T J-oo \ 2 c^(r) 



Wfi ~ (W(t)) 
1 f a w (r)V2 d£exp(-£ 2 

/7T J-oo 



(30) 



2 2 V(7 W (r)V2 



where, the error function is defined 



ll| as, 
2 



For large r, we have Callen-Welton theorem [6| which tells us 



er f( x )= f d£exp(-£ 2 ) (31) 
\/7T 7n 

16] 



Wr))-^ = l/3a 2 y (r). (32) 
Substituting for dissipation in terms of fluctuations in Eq. ( l30l) we get 

Note erf(oo) = 1 and erf (0) = 0. From Eq. ( 1331) we see that p(r) increases when r increases. 
In the limit of r — > oo, pfr) equals one-half. 

n 

This result can be understood as follows [7|. In the quasistatic limit of r — > oo, p(W; r) 
becomes more and more sharply peaked such that (W) — >■ Wr and cr^ — >■ 0. By Callen- 
Welton theorem Wj oc a w . This implies that aw 3> as —> 0. Thus and the 
aw both tend to zero, with aw 3> Wa- In other words — > faster than aw — > 0. This 
leads to p(r) — > \ in the reversible limit. 

However, dissipation is defined as Wd = (W) — Wr. In the reversible limit of r — > oo, we 
have (W) = Wr and hence dissipation is zero. 
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VI. WORK FLUCTUATIONS IN LIQUID CRYSTALLINE SYSTEM 



Let us now illustrate the above on a lattice model of a liquid crystalline system. To this 
end, we consider an L x L x L cubic lattice. Each lattice site holds a headless spin. The 
spins on nearest neighbor lattice sites interact with each other via Lebwohl-Lasher potential 



12j see below. Besides, each spin interacts with the external electric field, which is taken as 
the switching parameter. Without loss of generality we take the external field, to be in the 
z axis direction and switch its magnitude. The Hamiltonian for the nematic system is given 



by, 



H = -J]TP 2 (cos# M ) - ^]TP 2 (cos^). (34) 



In the above, the symbol (i, j) denotes that % and j are nearest neighbor lattice sites. The sum 
runs over all distinct pairs of nearest neighbor sites. We have employed periodic boundary 
conditions along x, y and z coordinates. 9ij denotes the angle between the spins and J 
measures the strength of interaction. In the simulation we express energy in units of J 
and hence set J = 1. E is the amplitude of the external electric field and 0$, the angle 
between the spin at lattice site i and the external electric field. P 2 (^) is the second Legendre 
polynomial given by 

PM = S -^f^± (35) 

We carry out Monte Carlo simulation of the response of the system to a process of switching 
E from say Eq to a value of Ef. 

VII. MONTE CARLO SIMULATION OF THE SWITCHING PROCESS 

Start with an arbitrary initial microstate (spin configuration) and employing Metropolis 



algorithm [13j, equilibrate the system at the desired value of (3 with E = E . Select randomly 
a microstate from the equilibrium ensemble. Calculate the energy. Keep the micro state 
the same and switch the field from E to E\ = E + AE. This is called the work step. The 
resulting change in energy, called work is denoted by Then implement a heat step via 
Metropolis algorithm over one Monte Carlo sweep as described below. 

Select randomly a spin and change its orientation by a random amount. This can done 



by the procedure suggested by Barker and Watts 



141 ] . Call the resultant microstate as 



trial state. If the energy of the trial state is less than that of the current state accept the 
trial state as the next microstate in the Markov chain. If the trial state energy is more, 
then calculate p = exp[— /3Ae] where Ae is the energy of the trial state minus that of the 
current state. Draw a random number (uniformly and independently distributed in the 
range zero to unity) and if it less than p accept the trial state; otherwise reject the trial 
state and take the current state as the next state in the Markov Chain. Carry out the above 
rejection/acceptance step L 3 number of times and this constitutes a Monte Carlo sweep. 
The heat step is followed by a work step on the micro state obtained at the end of the heat 
step. During the work step the field changes from E 1 to E 2 = E 1 + AE = E + 2AE. The 
work done is u 2 - The work and heat steps are repeated until E attains a pre - determined 
value, say Ef. 

Let n be the number of work steps required to switch the field from an initial value of E 
to a final value Ef. In other words E n = Ef = E + n AE, where AE = (Ef — E )/n. The 
switching time r is thus given by n. We have 

n 

W = I>. ( 36 ) 
i=i 

We carry out the simulation independently for a large number of times with the same 
switching protocol and construct a work ensemble {Wj} from which all the required statistics 
described below are calculated. In the simulation we have taken E = and Ef = 2. The 
results presented in this paper are for L = 3. Jarzynski identity is valid even for systems 
which are small. For a small system there would be large fluctuations and hence illustrating 
various issues, like relating work fluctuations to free energy differences and to dissipation 
becomes an easy task. This is the reason we have kept the system size small in this work. 
The size of the Monte Carlo ensemble of work values generated is one million. 

VIII. RESULTS AND DISCUSSIONS 

Fig. (1) depicts results on free energy difference (or reversible work) calculated from 
Jarzynski identity function of r with f3 — 1.5. 

We have also plotted in the same graph the free energy change from Callen-Welton theorem. 
As expected Callen-Welton theorem does not predict free energy change correctly for small 
switching time when the system is driven far from equilibrium. As r increases the results 
from Callen-Welton theorem converges to that from Jarzynski identity. 
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3 °500 500 1000 1500 2000 2500 3000 3500 
Switching time (x ) 

FIG. 1. Free energy difference AF versus switching time r. The upper curve (squares) denotes 
AF J calculated from Jarzynski identity. The lower curve (circle) denotes AF CW , calculated from 
Callen-Welton theorem. 

In Fig. (2) we have plotted dissipation given by Eq. ( TT8|) and denoted as IF/, as a function 
t. This has been calculated from the ensemble of work values generated by the Monte Carlo 
simulation of the switching process. The dissipation is large for small r. This is because the 
system is driven far from equilibrium. However as the switching time increases the process 
becomes more and more quasistatic. Dissipation decreases. In the limit of r — > oo it goes 
to zero. For large r dissipation can also be calculated from fluctuations employing Callen- 
Welton theorem. We have W$ w = /3a^/2. The inset in Fig. (2) shows W£ w - W J d . This 
quantity is large for small r. It decreases with increases of r and eventually goes to zero. 

We have depicted in Fig. (3) the probability of violation of the second law, calculated 
from the Monte Carlo ensemble of work values. Let p MC {r) denote this quantity which is 
calculated as follows. 

We count how many of the switching experiments give rise to a value of W less than Wr. 
This number divided by the total number of switching experiments carried out gives p MC {r). 
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FIG. 2. Dissipation defined as = (W) — Wr calculated from Jarzynski identity. The inset 
shows Wj W — W£ , where W^ w is dissipation calculated from C alien- Welton theorem. 

We find that p M (r) increases with increase of r as expected from the analytical results 
discussed in section V. 

The work distributions for representative values of r are depicted in Fig. (4). For small 
values of r we find the distribution is broad. As r increases the distribution becomes more 
and more sharply peaked. The distribution is expected to be Gaussian see 



15 



16]. 



IX. CONCLUSIONS 



Recent developments in non-equilibrium statistical mechanics embodied in various fluc- 
tuation theorems give us an insight into the foundational aspects of statistical mechanics 
and thermodynamics. There are several other methods which provide us a tool to calculate 
equilibrium quantities from non-equilibrium measurements. For example Sadhukhan and 
Bhattacharjee 17] have shown that Barkhausen noise process, repeated many times give ad- 
equate data to construct, in conjunction with work fluctuation theorem, a special matrix 
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FIG. 3. The probability of second law violation for various value of t . calculated from the Monte 
Carlo ensemble 



whose principal eigenvector provides equilibrium distribution. Transient entropy-fluctuation 
theorem of Evans and Searles jl9] . heat fluctuation theorems of Crooks[3l and the fluctu- 
ations theorems of Gallavotti-Cohen 



20| for steady state systems also can be employed to 
estimate equilibrium quantities from non-equilibrium measurements. In this paper we have 
employed work fluctuation theorem for estimating free energy differences. We have carried 
out Monte Carlo simulation of the lattice model of liquid crystalline system. Monte Carlo 
is best suited for lattice spin models. If the Hamiltonian involves continuous degrees of 
freedom then Molecular Dynamics simulations would be appropriate. We find that dissipa- 
tion defined as the excess of work done on the system over equilibrium free energy change 
tends to zero in the asymptotic limit of r — > oo. p(r) is given by the area under the work 
distribution in the tail region extending from — oo to Wr. We have shown analytically that 
the value of p increases with increase of r. In the quasistatic limit of r — >• oo, p{r) = 1/2. 
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FIG. 4. Distribution of work for r = 1, 16, 64, 256, 800 and 3200, from right to left. 
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